# misc.funLB -- 
#
# Copyright (c) 1988,2001,2003 Willem Jan Zaadnoordijk

from math  import exp,log,pi,atan
from funE_B import erfc


def ieer(x,a,b) :
	"function for transient line elements"
#  function returns the value of the integral of the e power of minus
#  x squared times the difference of the complementary error functions
#  of A*X and B*X evaluated from X=X to X=infinity
	sqrpi=1.77245385090551
	cutoffx=10.
	if (x > cutoffx) :
#     for large x the integral is zero because the exp(-x*x) part of
#     the integrand becomes zero and the rest stays finite
		rfieer=0.
	else :
#     the function can be expressed in terms of the funtion H
#     defined in eq.12 in Litkouhi and Beck (1982)
#     RFIEER(x,a,b)=.5/sqrt(pi)*[H(b*x,1/b)-H(a*x,1/a)]
#     without combining anything (which would speed things up in some
#     cases) this can be calculated as:
		rfieer=.5*sqrpi*( lb82(x,b)-lb82(x,a) )
	return rfieer


def lb82(rxin,rain) :
	"function for ieer"
#	RF..................L............B.......82
#	function taken from Litkouhi and Beck (1982)
#	equations (20) and (15), (16), (17a)
#	RFLB82(RXIN,RAIN)=H(RXIN*RAIN,1/RAIN)=H(RHX,RHP)=H(X,p)
	r2opi=2/pi
	rhpxlarge=10
	rhpxsmall=0.0001
	rhxbig=4
	rhp1plus=1.3
	rhp1minus=0.6
	rasmall=1.e-30
#	work with absolute values because of symmetry relations
	rx = abs(rxin)
	ra = abs(rain)
	rhx= ra*rx                                     # X of H(X,p)
	if ra < rasmall :
		rhp = 1/rasmall
	else :
		rhp = 1/ra                                   # p of H(X,p)
	rhpx= rx                                        # =pX

	if rhpx > rhpxlarge :
#		pX is large
#		this is never reached when called from RFIEER
#		where this condition is taken care of without
#		calling this routine RFLB82
		rflb82=0
	elif rhpx < rhpxsmall :
#		pX is small
		if rhp > 1 :
#			Litkouhi & Beck eq.21a: H(X,p)=2/pi*(arctan(1/p)-pX
			rflb82=r2opi*atan(ra)-rhpx
		else :
#			Litkouhi & Beck eq.21b:
#			H(X,p)=1-erf(px)erf(x)-2/pi*(arctan(1/p)-pX
			rflb82=1-(1-erfc(rhpx))*(1-erfc(rhx))  \
				-r2opi*atan(ra)-rhpx
	else :
#		pX is neither large nor small
		if rhx > rhxbig :
#			X is large
#			Litkouhi & Beck: equation 20
			rflb82=erfc(rhpx)
		else :
#			X is not large
			if rhp > rhp1plus :
#				p is decisively greater than 1
#				Litkouhi & Beck: equation 16
				rflb82=r2opi*exp(-rx*rx)*lbsu(rx,ra)
			elif rhp < rhp1minus :
#				p is decisively less than 1
#				Litkouhi & Beck: equation 15
				rflb82=1-(1-erfc(rhx))*(1-erfc(rx))- \
					r2opi*exp(-rhx*rhx)*lbsu(rhx,rhp)
			elif rhp == 1 :
#				p is exactly 1
#				Litkouhi & Beck: equation 17a
				rflb82=.5-.5*(1-erfc(rx))**2
			else :
#				p is about 1
#				Barnes & Strack ICAEM2003
				rflb82=barst(rhpx,rhx,rhp)

#	H(X,-p) =-H(X,p)
#	H(-X,p) = H(X,p)
#	H(a(-x),1/a) = H(-ax,1/a) = H(ax,1/a)
#	H(-ax,1/(-a))= H(-ax,-1/a)= H(ax,-1/a) = H(ax,1/a)
	if rain < 0:
		rflb82=-rflb82
	return rflb82


def lbsu(rx,ra) :
	"function for lb82"
#	                      (-1)**n * e_n(RX**2) * RA**(2n+1)
#	sum for n=0 to inf of ---------------------------------
#	                                    2n+1
#	used for formula (15) and (16) of Litkouhi and Beck (1982)
#	e_n from Abramowitz and Stegun: page 262, equation 6.5.11
	nordmx=100
	racc=1.e-8
	rcut=50
	if rx > rcut :
		rflbsu = 0
	else :
		rsign = 1
		rxnonf= 1
		renrx = 1
		raa   = ra**2
		ra2np1= ra
		r2np1 = 1
		rsum  = ra
		r0    = abs(ra)*racc
		i     = 1
		while i <= nordmx :
			rsign = -rsign
			rxnonf= rxnonf*rx/i
			renrx = renrx+rxnonf
			r2np1 = r2np1+2
			ra2np1= ra2np1*raa
			ri    = rsign*renrx*ra2np1/r2np1
			rsum  = rsum+ri
			if abs(ri) >= r0 :
				break
			i = i+1
		rflbsu=rsum

	return rflbsu


def barst(rhpx, rhx, rhp) :
	"function for lb82"
#	approximation for H(X,p) function in the range
#	0.0001 <= pX <= 4, X<=4, and 0.6 <= p <= 1.3
#	the function H(X,p) is defined by Litkouhi and Beck (1982)
#	the approximation was formulated by Barnes & Strack in
#	a paper for the ICAEM conference in France in April 2003
#	which was cancelled unfortunately
	r2opi=2/pi
	nmax=100
	rtol=1.e-8

	rtest  = brf(rhx)
	rphia  = 0.5*rtest
	rphib  = rtest**2
	rsum   = rphib
	rtest  = abs(rphib)*rtol
	rhxx   = rhx**2
	rhpp   = rhp**2
	r1mppo2=.5-.5*rhpp
	r1pppn =1
	i = 1
	while i <= nmax :
		rphic = ( (2*i-1-2*rhxx)*rphib   \
			         + 4*rhxx*rphia - 0.5 \
			     ) / i
		r1pppn= r1pppn*r1mppo2
		ri    = r1pppn*rphic
		rsum  = rsum+ri
		if abs(ri) < rtest :
			break
		rphia = rphib
		rphib = rphic
		i = i+1
	rfbarst = erfc(rhpx)-r2opi*rhp*exp((-1-rhpp)*rhxx)*rsum
	return rfbarst


def brf(rhx) :
	"function for barst"
#	function from Barkley Rosser (1948): equation 6-2
#	that is used in the approximation for H(X,p) by
#	Barnes & Strack in a paper for the ICAEM conference
#	in France in April 2003 which was cancelled unfortunately
	rsqrpi=1.772453850905516027298167
	rsqrpio2=rsqrpi/2
	rfbrf = exp(rhx**2)*rsqrpio2*erfc(rhx)
	return rfbrf
